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ABSTRACT 

The cosmic infrared background (CIB) includes roughly half of the energy radiated by all galaxies at all wavelengths across cosmic 
time, as observed at the present epoch. The PACS Evolutionary Probe (PEP) survey is exploited here to study the CIB and its redshift 
differential, at 70, 100 and 160 /jm, where the background peaks. Combining PACS observations of the GOODS-S, GOODS-N, 
Lockman Hole and COSMOS areas, we define number counts spanning over more than two orders of magnitude in flux: from ~ 1 
mJy to few hundreds mjy. Stacking of 24 yum sources and P{D) statistics extend the analysis down to ~0.2 mjy. Taking advantage 
of the wealth of ancillary data in PEP fields, differential number counts d 2 N/dS jdz and CIB are studied up to z = 5. Based on 
these counts, we discuss the effects of confusion on PACS blank field observations and provide confusion limits for the three bands 
considered. While most of the available backward evolution models predict the total PACS number counts with reasonable success, the 
consistency to redshift distributions and CIB derivatives can still be significantly improved. The new high-quality PEP data highlight 
the need to include redshift-dependent constraints in future modeling. The total CIB surface brightness emitted above PEP 3cr flux 
limits is vl,, = 4.52 ± 1.18, 8.35 ± 0.95 and 9.49 ± 0.59 [nW nr 2 sr 1 ] at 70, 100, and 160 yum, respectively. These values correspond 
to 58 ± 7% and 74 ± 5% of the COBE/DIRBE CIB direct measurements at 100 and 160 yum. Employing the P(D) analysis, these 
fractions increase to ~65% and ~89%. More than half of the resolved CIB was emitted at redshift z < 1. The 50%-light redshifts lie 
at z = 0.58, 0.67 and 0.73 at the three PACS wavelengths. The distribution moves towards earlier epochs at longer wavelengths: while 
the 70 /im CIB is mainly produced by z < 1.0 objects, the contribution of z > 1.0 sources reaches 50% at 160 yum. Most of the CIB 
resolved in the three PACS bands was emitted by galaxies with infrared luminosities in the range 10" - 10 12 L G . 

Key words. Infrared: diffuse background - Infrared: galaxies - Cosmology: cosmic background radiation - Galaxies: statistics - 
Galaxies: evolution 
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' 1. Introduction sources, such as active galactic nuclei, are expect ed to be only 



H | 5-20% to the total EBL in the mid- and far-IR ( e.g. Matut e et all 

TO , With the exception of the cosmic microwave background, which 2 006; Jau zac et alJl2()TH iDraner & Ballantvnell201 lb . The onti- 

represents the relic of the Big Bang, emitted at the last scatter- cal and i n f ra red backgrounds dominate the EBL by several or- 

ing surface, the extragalactic background light (EBL) is the inte- ders f magnitude with respect to all other spectral domains, 
gral of the energy radiated by all galaxies, from y-rays to radio 

frequencies, across all cosmic epochs. Its energy density distri- The cosrmc infrared background (CIB) was detected and 

bution is characterized by two primary peaks: the first at X ~ 1 measured for the first time in the mid nineties, analyzing the 

fim produced by starlight, and a second peak at ~ 100 - 200 data obtained with the DIRBE and FIRAS instruments aboard 



fim mainly du e to starlight absorbed and reprocessed by dust in the C o smic Background Explorer (COBE) satel li te dPueet et al 



galaxies (e.g. lHauser & Dwek 2001). Contribution from other 



19961: lHauser et alJ Il998t iFixsen et alJ Q.998; Lagache et al 



1999i 120001) . In the global budget, the CIB provides roughly 



Send offprint requests to: Stefano Berta, e-mail: berta@mpe.mpg . de half o f the total EBL (e.g. |Dole et al.| |2006|; |Hauser & Dwek[ 

* Herschel is an ESA space observatory with science instruments pro- |200l|). Since in the local Universe the infrared output of galax- 

vided by European-led Principal Investigator consortia and with impor- ies is only a third of the em ission at optical wavelengths (e.g. 

tant participation from NASA. ISoifer & Neugeba uedll99ll) . this implies a strong evolution of 



1 



S. Berta et al.: Building the cosmic infrared background brick by brick with Herschel/PEP. 



infrared galaxy populations, towards an enhanced far-IR output 
in the past, in order to account for the total measured CIB. 

The discovery of large numbers of distant sourc es emitting 
a substantial amount of their energy in the IR (e.g. | Smail et all 
ll997tlAussel et al.ll999LlElbaz et al.ll999tlPapovich et all2 004. 
among many others) demonstrated that, although locally rare, 
powerful IR galaxies are indeed numerous at high redshift. The 
deep extragalactic campaigns carried out in the nineties and 
2000's with the Infrared Space Observatory (ISO, see Genzel & 
Cesarsky 2000 for a summary), and the Spitzer Space Telescope 
(Werner et al. l2004t e.g. Papovich et al. |2004) were very effi- 
cient in identifying and characterizing large samples of mid-IR 
sources. In contrast, at wavelengths near the CIB peak (100-200 
pm), the performance of these telescopes was strongly limited 
by their small apertures (< 85 cm diameter), the prohibitive con- 
fusion limits, and detectors' sensitivity. Spitzer surveys at 70 and 
160 pm produced limited samples of distant far-IR objects (e.g. 
iFraver et al.ll2009h : in the 160 pm Spitzer/MIPS band only ~ 7% 
of the CIB w as resolved into individually detected objects (Dole 
et al., 2004, see Lagache et al. 2005 for a review). Performing 
stacking of 24 pm sources. lDole et al] (12006) retri eved more than 
70% of the far-IR background at 70 and 160 pm. Bethermin et al. 
(1201 Oal) extended this analysis by extrapolating to very faint flux 
densities using a power-law fit to stacked number counts, and 
produced an estimate of the CIB surface brightness in agreement 
with direct measurements. 

Launched in May 2009, Herschel dPiibratt et al.l2010h is pro- 
viding stunning results: its large 3.5 m mirror, and the high sensi- 
tivity Photodetector Array Camera & Spectrometer (PACS, per- 
forming imaging at 70, 100, 160 pm; Poglitsch et al.. 1201 Oh are 
tailored to overtake the confusion and blending of sources that 
were hampering the detection of fain t far-IR sources in previous 
space missions. In lBerta eta D (120101 hereafter called "Paper I"), 
we exploited data from the PACS Evolutionary Probe (PEP) 
survey, covering the GOODS-N, Lockman Hole and COSMOS 
(to partial depth only) fields, and we resolved into individual 
sources 45% and 52% of the CIB at 100 and 160 pm. At longer 
wavelengths, in the spectral domain cove red by the Spectra l 
and Photometric Im aging Receiver (SPIRE. IGriffin et alJ uOlO). 
lOliver et all d2010h resolved 15%, 10% and 6% of the CIB at 
250, 350 and 500 pm. These fractions incr eased to 64%, 60% 
and 43 %, respectively, us i ng a P( D) analysis dGlenn et alJ2010h . 
Finally, iBethermin et al.l d2010bl) retrieved roughly 50% of the 
CIB applying stacking of 24 pm sources to BLAST 250, 350 
and 500 pm maps. 

The PACS Evolutionary Probe (PEfQ) is one of the ma- 
jor Herschel Guaranteed Time (GT) extragalactic projects. It is 
structured as a "wedding cake" survey, based on four different 
layers in order to cover wide shallow areas and deep pencil- 
beam fields. PEP includes the most popular and widely studied 
extragalactic blank fields: COSMOS (2 deg 2 ), Lockman Hole, 
EGS and ECDFS (450-700 arcmin 2 ), GOODS-N and GOODS- 



ing detail, from its initial discovery with COBE to the Spitzer 
era. The real conundrum is now represented by how the CIB 
evolves as a function of cosmic time. This information is not 
only important to constrain models of galaxy evolution, but also 
to shed light on the intrinsic spectra of TeV sources, such as 
BLAZARs, whose y-ray photons interact with the CIB (e.g. 
Dommguez et al. 201 1; Franc eschini et al.ll2008tlMazin & Raud 
2007). Thanks to Herschel capabilities, in Paper I we have ini- 
tiated the study of the redshift build-up of the local far-IR CIB, 
by associating to each object detected by PACS in GOODS-N 
a complete UV-to-FIR SED and splitting number counts and 
CIB into redshift slices. Hauzac et al] d201 lb performed a similar 
study based on stacking on 70 and 160 pm Spitzer/MIPS maps of 
COSMOS, and lLe Floc'h et al.l(l2009l) performed a similar anal- 
ysis at 24 pm in the same field. Here we take advantage of deeper 
PACS observations in GOODS-S and full coverage in COSMOS 
to further explore the history of CIB, with the aim to understand 
when it was primarily emitted. 

In this paper, the cosmic far-IR background is reconstructed 
"brick by brick" in three main "storeys". In Section [2] the re- 
solved component of number counts is presented, including slic- 
ing in redshift and a comparison to available evolutionary mod- 
els. Section [3] deals with stacking of 24 pm sources onto PACS 
maps. The third step (Sect.|4]i extends the analysis to even fainter 
fluxes using the "probability of deflection", P(D), technique. 
After a brief digression about source confusion in Sect. [5] we 
discuss the observed properties of CIB in Sect. [6] Finally, Sect. 
[7]draws a summary of our findings. 

2. Level 1 : detected sources 

In Paper I we studied number counts based on Science 
Demonstration Phase (SDP) and early routine phase data, i.e. 
based on the GOODS-N, Lockman Hole and COSMOS fields. 
At the time of SDP, COSMOS wa s available only to 2/3 of its 
nominal depth. lAltieri et all (12010) exploited gravitational lens- 
ing in the Abell 2218 cluster of galaxies to estimate the lens 
amplification of the fluxes of background galaxies. In this way, 
they were able to push the analysis down to 1 mJy at 100 and 
160 pm, thus breaking the confusion limit for blank field obser- 
vations (see Sect. [5J. 

Here we complete the Herschel/PACS view of far-IR number 
counts and CIB, making use of the deepest PEP field, GOODS- 
S, and the full coverage of COSMOS, in addition to Paper I 
results. Data reduction, catalogs construction and simulations, 
aimed at deriving completeness, fraction of spurious sources and 
photometri c reliability, are described in Lutz et al. (in prep.) and 
Bert a et al. | d2010b . Table Q] summarizes the main properties of 
the fields taken into account here. 

We applied the method described bv lCharv et al.l d2004l) and 
Smai l et al.l d!99 5) to correct number counts for incompleteness, 
on the basis of simulations. The distribution of input and out- 



S (-200 arcmin 2 ). In addition, the fourth tier of the "cake" con- put fluxes in simulations is organized in a matrix P u , so that i 



sists of ten nearby lensing clusters, offering the chance to break 
the PACS confusion limit thanks to gravitational lensing. An 
in depth description of PEP fields and survey properties is pre- 
sented by Lutz et al. d201 ll in prep.). 

Here we extend the analysis carried out in Paper I to the 
deepest field observed by PEP, GOODS-S, reaching a 3 <x 
depth of 1.2 mJy at 100 pm, and including also the 70 pm 
PACS bandpass. The surface brightness of the CIB and its spec- 
tral energy distribution (SED) has become known with increas- 



http://www.mpe.mpg.de/ir/Research/PEP/ 



represents the input flux and j the output flux. In other words, 
the zj-th element of the matrix gives the number of sources with 
z'-th input flux and j-th output flux. The way is built implies 
that TjjPij ^ 1 represents the completeness correction factor 
at the z'-th input flux in simulations. In order to correct the ob- 
served counts, Pij is re-normalized such that the Py equals 
the number of real sources detected in the j-th flux bin. Finally, 
the completeness-corrected counts in the given z'-th bin are given 

by ^ p renorm 

Figure Q] presents PACS number counts at 70, 100 and 
160 pm, normalized to the Euclidean slope (i.e. the slope ex- 
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Fig. 1. Differential number counts in the three PACS bands, normalized to the Euclid ean slope (dN/dS oc S 2 S ). Filled/open 
symbo l s belong to flux bins abov e /below the 80% comple te ness limit. Mo d els be l ong to [ Lagache et alj d2004l). Franceschini et al.l 
(120101). iRowan -Robinsonl d2009h . EeBorgne et"al] d2009t) . IVahante et all d2009h . lLacev et all d2010t) . [BetherminltalJ (|2010cJ), 
Mars den et al.l d2010|). Gruppioni et al. d201 ll in prep.), Niemi et al. (1201 jl in pre p.1. Shaded areas represent ISO and Spitzer data 
(Rodighiero & Franceschini 2004; Heraudeau et al. 2004; Bethermin et al. 2010a); hatched areas belong to Spitzer 24 fim stacking 
(Bet hermin et alJuOlOal) . Left and right panels present individual fields and averaged counts, respectively. 



pected for a uniform distribution of galaxies in Euclidean space, 
dN/dS oc S~ 2 5 ). Error bars include Poisson statistics, flux cal- 
ibration uncertainties, and photometric errors. The latter have 
been propagated into number counts via 10 4 realizations of ran- 
dom Gaussian flux errors applied to each PACS source, using a 
dispersion equal to the local measured noise. It is worth to note 



that in most cases the faint end of counts derived in shallow fields 
is consistent with results from deeper fields, thus confirming the 
validity of completeness corrections. Number counts, their cor- 
responding uncertainties and completeness values are reported 
in Tabs. [3] S and|5] 
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Fig. 2. Field-to-field variations of number counts in COSMOS. Top panels: number counts split into sub-areas (lines) compared 
to full-field counts (black dots). Bottom diagram: standard deviation among the sub-fields (top), as a function of flux. In the very 
bottom panels, we show the comparison between the field-to-field cr(var) and the uncertainty o> obtained combining Poisson 
statistics, photometric errors and systematics. Shaded areas represent lcr uncertainties on o~(var) and cr(var)/crp determinations. 



The left panels in Fig. Q] show the counts for each field sep- 
arately. Filled symbols represent data above 80% completeness, 
while open symbols extend down to the 3<x detection threshold. 
The grey shaded areas in Fig. [1] belong to pre-Hersche l num- 
ber counts. We include the [kodighiero & Franceschini (2004) 
and iHeraudeau etail ((2004) 95 pm IS O data, and th e 70 pm 
and 160 pm Spitzer number counts by Bethermin et al. (2010a), 
based on GOODS/FIDEL, COSMOS and SWIRE fields. Both 
individual detection and stacked Spitzer counts (hatched areas) 
are shown. PACS counts and previous results are in good agree- 
ment, over the flux range in common. The PEP deepest field, 
GOODS-S, extends the knowledge on far-IR number counts one 
order of magnitude deeper in flux than Spitzer individual de- 
tections at 160 pm and roughly 5 times deeper at 70 pm. The 
3cr limit in GOODS-S (1.2 mJy and 2.4 mJy at 100 and 160 
Ai m, respectively) is v ery close to the effective depth reached 
bv lAltieri et alJ d2010h in the Abell 2218 cluster when studying 
lensed background galaxies, but our improved source statistics 
provide much tighter uncertainties. 



In order to provide a single reference counts description, the 
number counts belonging to the four PEP fields studied have 
been combined via a simple average, weighted by their respec- 
tive uncertainties in each flux bin. Results are included in Tabs. 



[3] |H [3J and are shown in the right panels of Fig. [T] as compared 
to a collection of model predictions (see Sect. 12.41 ). 

The depth reached by PACS/PEP allows us to accurately 
probe the faint-end of counts. The peak in the normalized counts 
is well sampled and turns out to lie at ~4 mJy at 70 pm, ~ 10 mJy 
at 100 pm, and ~20-30 mJy at 160 pm. The differential counts 
are reproduced by a broken power law (dN/dS oc S a ), charac- 
terized by a break at flux S break and two distinct slopes at the at 
faint/bright sides of the break. A weighted least squares fit was 
performed on the data, and the results are presented in Tab.|2l for 
different fields and flux ranges. Breaks happen at ~3.5 mJy at 70 
pm, ~5.0 mJy at 100 pm, and ~9.0 at 160 pm. Uncertainties at 
the bright end are dominated by Poisson statistics, and nearly- 
Euclidean slopes are allowed. 

2. 1. Field to field variations 

The large area probed by COSMOS (~2 deg 2 ) can be used to 
test the effect of field-to-field density variations on the bright- 
end of number counts. We adopt a fully empirical method, de- 
scribing the variance in number counts coming from the inferred 
variations, while a full clustering quantification goes beyond the 
scope of this paper. To this aim, we split the COSMOS field into 
a number of fully independent tiles, probing different angular 
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scales: 28 tiles with size -200 arcmin 2 , similar to the size of 
GOODS fields, 12 Lockman Hole like tiles (~500 arcmin 2 ), and 
4 tiles of ~1500 arcmin 2 each. 
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Fig. 3. Comparison between photometric and available spectro- 
scopic redshifts in GOODS-S, GOODS-N and COSMOS, re- 
gardless of PACS detections. Dotted and dashed lines mark 10% 
and 20% uncertainty levels. 



Number counts in each tile were computed as described be- 
fore and are plotted in Fig. [2] as thin solid lines. Black symbols 
represent the average counts in each flux bin, and the associ- 
ated error bars are computed in the usual manner accounting for 
Poisson statistics, systematics and photometric uncertainties. 

The properties of field-to-field standard deviation are stud- 
ied in the bottom diagrams of Fig. [2] As expected, the <j(var) 
in number counts increases as a function of flux (upper panels), 
both because more luminous sources are rarer in the sky and 
because lower-redshift objects are dominating at these fluxes, 
hence probing a smaller volume. Over the whole flux range 
covered by this analysis, the field-to-field deviation is compara- 
ble to the uncertainty crp obtained combining Poisson statistics, 
photometric errors and systematics (bottom panels). Overall, 
<r(var)/o-p is slightly larger than unity within the errors. This ef- 
fect can be explained considering that neighboring sub-tiles are 
not fully independent: because of clustering, a non-negligible 
correlation term contributes to cr(var). At larger scales (dotted 
green lines), <r(var)/o-p is more noisy because of the limited 
number of tiles available. 



2.2. Ancillary data and multi-wavelength catalogs 

The fields observed with PACS as part of the PEP survey ben- 
efit from a plethora of ancillary data, spanning from the x-rays 
to radio frequencies. We took advantage of these data to build 
reliable multi-wavelength catalogs and associate a full spectral 
energy distribution (SED) and a redshift estimate to each PACS- 
detected object. 

As described in Paper I, a PSF-matched catalog was created 
in GOODS-N, including photometry from GALEX far-UV to 
Spitzer IRAC and MIPS 24 fim. The Southern GOODS field is 
rich in covera ge as well. Here we adopt the PSF-matched cat- 
alog built by Grazi an et alJ d2006h. t o which we add the 24 //m 
photometry by Mag nelli et alJ (120091) and a collection of spectro- 
scopic redshif t s for more than 3000 source s (Balestra et al. 2010; 
Pope sso et all 120091: [Santini et al] l2009t IVanzeila et all 2008: 
Le Fevre etap 120051: [Mignoli et al.l 120051: iDohertv et al] 1 2005 
Iszokolv et al] l2004t jPickinson et al] ' 



2004; van der Wei et al 



2004; Stanwav et al. 2004b a; Stroker et al. 2004; Bunker et al. 
120031: ICroom etal.1 1200 U ICristiani et all 120001) . Finally, we 
browsed the COSMOS public databas^ and com bined the U- 
to-K broad- and inter mediate-band photometry (Capa k~et al.l 
120071: lllbert et al]|2009l containing 2 ,017,800 sources), the pub- 
lic IRAC catalogs, the 24 urn data dLe Floc'h et al]|2009l) and 
the available pho t ometric (Ilbert et al. 2009) and spectroscopic 
dLillv et al.l 120091: iTrump et al] 120091) redshifts. As far as the 
Lockman Hole is concerned, an extensive ancillary catalog is 
currently on the make by Fotopoulou et al. (1201 ll in prep.), but 
is not available at the time of this analysis, hence this field will 
not be used in this piece of analysis. 

PACS catalogs were matched to the ancillary source lists 
by means of a multi-band maximum likelihood procedure 
dSutherland & Saunderslfl992l) . starting from the longest wave- 
length available (160 fim, PACS) and progressively matching 
100 //m (PACS), 70 fim (PACS, GOODS-S only) and 24 jum 
(Spitzer/MIPS) data. 

When no spectroscopic redshift is available, a photometric 
estimate is necessary. In GO ODS-S we make use of the avail- 
able photometric redshifts by iGrazian et al .1 d2006l) . while new 
photo-z's were produced in GOODS-N, exploiting the wealth 
of multi-wavelength data collected, and adopting the EAZY 
(Bra mmer et al.ll2008l) code. Up to 14 photometric bands were 
used, depending on the data available. The top panel in Fig. [3] 
presents the comparison between photometric and the available 
spectroscopic redshifts. The fraction of outliers, defined as ob- 
jects having A(z)/(1 + z spe c) ^ 0.2, is ~6% over the whole sam- 
ple of spectroscopic redshifts, and decreases to ~2% for sources 
with a PACS detection. Most of these outliers are sources with 
few photometric points available, or SEDs hardly reproduced by 
the available templates. The median absolute deviation (MAE0) 
of the A(z) distributiorflis 0.040 for the whole catalog, and 0.038 
for PACS-detected objects with spec-z available. 

In the area covered by ancillary data, roughly 60-65% of 
GOODS-N sources detected by PACS in either band have a 
spectroscopic redshift estimate. In the GOODS-S MUSIC area 
(IGrazian e t al. 2006) this fraction is as high as ~80%. Roughly 
95% of these spectroscopic redshifts in GOODS-S lie at z < 2.0, 
with an almost complete coverage. For the remaining PACS 
sources we adopt photometric redshifts, obtaining a 100% red- 
shift completeness above the 80% photometric completeness 



2 http://irsa.ipac.caltech.edu/data/COSMOS/ 

3 MAD(x) = median (|x - median (x)\) 

4 where here A denotes the difference between photometric and spec- 
troscopic redshift. 
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malized to a 1 deg 2 area, and above the 80% photometric completeness limit. Black lines refer to models and are reported only at 
the GOODS-S and COSMOS depths for clarity sake. 



threshold (see Tab.Q]). As far as COSM OS is conce rned, the pub- 
lic photometric redshift catalog filbert et al.ll2009l) is limited to 
a magnitude / < 25, thus producing a redshift incompleteness in 
the PACS-selected sample. Only ~7 5% of PACS sourc es have a 
photometric redshift estimate in the lllbert et alj d2009l) catalog. 
This incompleteness is independent of PACS fluxes. We derived 
new photometric redshifts using the EAZY code and exploiting 
the public COSMOS datasets. The bottom panel of Fig.[3]shows 
the results. The fraction of outliers is ~1% in thi s case, and 
MAD(A(z)) =a 0.01. This result is similar to that bv lllbert etafl 
(2009) for the objects in common. 

Given the non-null fractions of outliers in Fig. [3] it is possible 
that the number of sources at high redshift (e.g. z > 2) in PEP 
catalogs is in part contaminated by "catastrophic photometric - 
redshift failures", mainly represented by low-z objects with 
wrong, high photo-z. Because of the paucity of z > 2 PACS 
sources benefiting from a spectroscopic follow-up, the only vi- 
able approach to test this effect is to assume that the distribution 
of A(z)/(1 + Zspec) of PACS-detected objects is similar to that of 
the general galaxy population. For this reason, the resulting con- 
tamination fraction has to be considered as an upper limit only. 

The fraction of potential contaminants at z > 2 is computed 
in two steps. First the fraction of z < 2 objects having z sp ec < 2, 
but A(z)/(1 + Zspec) ^ 0.2 and z p hot > 2 is derived; then this 
is re-scaled to the ratio of z % 2 PACS sources. It turns out 
that up to -25% of GOODS-S z > 2 PACS objects might have 
an improperly attributed high redshift. Similar results were ob- 
tained in GOODS-N, while this effect cannot be properly tested 
in COSMOS, because public z > 2 spectroscopic redshifts are 
currently still lacking. 

We also note that the opposite phenomenon — namely high- 
redshift sources with wrong photo-z potentially being shifted at 
low-z — does not play a significant role here, because the vast 
majority of z < 2 PEP/PACS objects benefits from a spectro- 
scopic redshift measurement. 



2.3. Contribution to the counts from different epochs 

We exploit the rich ancillary information described in the pre- 
vious Section to perform a detailed study of counts across cos- 
mic time. The redshift derivative dN/dz of PACS sources above 
the 80% photometric completeness limits, normalized to 1 deg 2 , 
is shown in Fig. [4] for the three fields and bands considered 
here. The covered flux ranges are quoted on each panel. In the 
COSMOS area, we sample the bright end of PACS counts. The 
distribution in this field peaks at z < 0.5 in both 100 and 160 /im 
bands, with 70-80% of all sources lying below z — 1 and -20% 
between z = 1 and z = 2. For the deeper GOODS-N data, the 
peak of the redshift distribution shifts to z ~ 0.7. Although the 
small sampled area limits source statistics, objects at high red- 
shift (up to z — 5) start to pop up. Our deepest field, GOODS-S, 
covered in all three PACS bands, displays some remarkable fea- 
tures. Overall dN/dz is now peaked around z ~ 1. It is possible 
to re cognize two wel l known structures at z — 0.7 and z - 1 . 1 
(e.g. iGilli et al.l l200l IVanzella et al.ll2005h . which produce nar- 
row and intense spikes in the distribution at 100 and 160 //m. 
On the other hand, at 70 /-im these structures are barely seen. At 
higher redshift, a broad "bump" is detected, between z = 2 - 3. 
This peak cannot be identified in the shallower 70 fim data, but 
is outstanding at the other two PACS wavelengths. Cutting the 
GOODS-S catalog at the depth reached by GOODS-N (5.5 mJy 
and 11.0 mJy at 100 and 160 fim, respectively), the high-z fea- 
ture disappears and the redshift distribution resembles that of 
GOODS-N, with only a few sources left above z > 2. Similarly, 
when cutting GOODS-S at the COSMOS 80% depth, we re- 
trieve a distribution peaked at z = - 0.5, obviously with much 
poorer statistics than in the COSMOS field itself. An extensive 
analysis of PACS GOODS-S large scale structure at z = 2 - 3 
and of a z — 2.2 filamentary over density is being presented by 
Magliocchetti et al. d20lU sub.). Table |6]reports dN/dz [deg -2 ], 
as derived in these three PEP fields. The median redshifts of the 
sources detected in GOODS-S are z ~ 0.6, ~ 0.7, ~ 0.8 in the 
three PACS bands, but would shift if the known spikes in dN/dz 
were not present. 
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Finally, it is possible to split number counts in GOODS-S, 
GOODS-N and COSMOS into redshift bins, similarly to what 
was done in Paper I for GOODS-N. Counts are then constructed 
in four broad bins, in order to allow for a sufficient number 
of sources in the small fields. In every redshift bin we apply 
the same completeness correction derived for the total number 
counts at the given flux. Similarly to what was done for the total 
counts, we build combined number counts through a weighted 
mean in each flux bin. The average number counts, sliced in red- 
shift intervals are reported in Tabs. [7] [8] and|9] and are shown in 
Fig-0 PACS observations in the three fields nicely complement 
each other: the southern GOODS field reaches deep flux den- 
sities, not probed previously, and also traces the high redshift 
population that was missing in GOODS-N. The bright end of 
counts is sampled by observations in the COSMOS area, which 
not only contribute to the low-redshift counts, but also nicely 
trace the bright component up to z — 2 and beyond. 
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2.4. Comparison to model predictions for counts across time 

Number counts and redshift distributions encode the evolution 
of galaxies: the upturn at intermediate fluxes and the over- 
Euclidean slope at the bright side of the peak are usually inter- 
preted as signatures of strong evolution in the properties of the 
underlying galaxy populations, and have stimulated a plethora 
of model interpretations. 

We are comparing here our results to examples of two 
basic classes of models attempting to reproduce these observ- 
ables. Backward evolutionary models transform the statistical 
properties of far-IR galaxies observed at known redshift (e.g. 
local luminosity functions) into observables at any redshift 
(e.g. galaxy counts, CIB brightness, etc.) assuming a library of 
template SEDs and parametric laws of luminosity and/or density 
evolution. These models do not implement fundamental physi- 
cal information but simp ly attempt to describe the evolution of 



2004; Francesc 



LeBorgne et al 



uni et al] 
2009t 



galax y p opulations (e.g.|L agache et al 
20 10t iRowan- Robinso nl 120091: 

Valiante et all l2009t iBethermin et al.1 120 lOd iMarsden et all 
2010; iGruppioni. C. Pozzi. E. Zamorani. G.. Vignali. C. et al.l 
201% . 

Very simple in their principles, backward models are 
thus parameterizations embedding the knowledge obtained 
from previous observations. The different adopted recipes 
were optimized to reproduce available observables such as 
mid-IR ISO and Spitzer number counts and far-IR Spitzer 
counts based on detections and — in some cases — 
stacking, sub-mm counts, as well as redshift distributions. 
A couple of models took advantage of early Herschel 
results, namely SDP PACS and SPIRE number counts 
(Beth ermin et al. 201 Pel) and PACS luminosity functions up t o 



z ~ 3 dGruppioni. C. Pozzi. R. Zamorani. G.. Vignali. C. et al] 
1201 ll) . Generally speaking, confronting backward models to the 
new, detailed Herschel data is a test to their flexibility in produc- 
ing reliable predictions at far-IR wavelengths. 

At the opposite side, forward evolution models simulate the 
physics of galaxy formation and evolution forward in time, from 
the Big Bang to present days. Current implementations are based 
on semi-analytic recipes (SAM) to describe the dissipative and 
non-dissipative processes influencing galaxy evoluti on, framed 
into A -CDM dark matter numerical simulations (e.g. Lace v et al] 
120101 Niemi et al. in prep.). Galaxy radiation, including far-IR 
emission, is computed with spectrophotometric synthesis and ra- 
diation transfer dust reprocessing, given the fundamental prop- 
erties of the galaxies in the model. 
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Fig. 5. Differential number counts in the three PACS bands, nor- 
malized to the Euclidean slope and split in redshift bins. The 70 
fim counts belong to the GOODS-S field, while those at 100 and 
160 fim were obtained via a weighted average between GOODS- 
S, GOODS-N and COSMOS. See Fig. Q] for models references. 



Figures Q] and [5] show a collection of models overlaid to the 
observed PACS counts. Most of the models reproduce fairly well 
the observed 100 and 160 fim normalized counts, on average the 
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most s u ccessful being Franceschini et al. ( 20101). Marsden et al.1 
(1201 Oh. iRowan-Robinsonl (l2009t). IValiante et al.1 (120091) and 
Gruppioni, C, Pozzi, F., Zamora ni. G.. Vignali. C. et alj 
(120111) . Very different assumptions produce relatively similar 
number counts predictions. The iFranceschini et al] {2010) 
model adopts four different galaxy populations, including 
normal galaxies, luminous infrared galaxies, AGNs and a class 
of strongly-evolving ULIRGs, which dominate above z > 1.5, 
but is negligibl e at later epochs, r esemb ling high-z sub-mm 
galaxies. Also Rowan-Rob insonl d2009l) uses four galaxy 
populations (cirrus-dominated quiescent galaxies, M82-like 
starbursts, Arp220-like extreme starbursts, AGN dust torii), but 
employs analy t ic ev olutionary functions without discontinuities. 
Valia nte et al.l (120091) describe galaxy population taking into 
account the observed local dispersion in dust temperature, the 
local observed distribution of AGN contribution to Ljir as a 
function of luminosity, and an evolution of the local lum i nosity- 
temperature relatio n for IR ga l axies. [Marsden et al. (2010) 
base their SEDs on iDraine & Lil (120071) prescriptions and tune 
them to reproduce the local color-dependent far-IR luminosity 
functions; they include luminosity, density and color evolutions 
to fit sub-mm counts, redshift distributions and EBL. Finally , 
iGruppioni. C.. Pozzi. E. Zamorani. G.. Vignali. C.. et alj d201 lb 
include a significant Seyfert-2 population based on a fit to 
Herschel LFs (see also lGruppioni ^taT]l20 1 Oh . 

Semi-analytical approaches surely represent a more com- 
plete view of galaxy evolution, including a wide variety of 
physics in a single coherent model. They cover a wide range 
of observational data: UV, optical, near-IR luminosity functions, 
galaxy sizes, metallicity, etc. Moreover, at far-IR wavelengths, 
even in the case that global properties such as star formation 
and AGN activity are correctly modeled, a further complication 
arises from the assumptions about dust content and structure, 
which need to be invoked. As a result, the large number of pa- 
rameters involved goes at the expense of inference precision: the 
performance of SAM models with respect to PACS observables 
needs still substantial tuning. 

As described above, it seems that pre-Herschel models, with 
only a few exceptions, are quite successful at reproducing total 
number counts despite the range and diversity of the employed 
solutions, thus showing that the discriminatory power of this ob- 
servable is rather limited. 

The right answer comes from the redshift information avail- 
able in the selected PEP fields. Figure [4] presents the redshift 
distribution dN/dz of PACS galaxies, and Fig. [5] shows number 
counts split in redshift slices (d 2 N/dS /dz). The combination of 
the two is a real discriminant: model predictions are now dra- 
matically different. 

None of the available models provides a convincing coherent 
prediction of the whole new set of observables (both dN/dz and 
d 2 N/dS /dz), over the flux range covered. One common source 
of discrepancy seems to be the well-know degeneracy between 
dust temperatu re (and hence lu minosity and total dust mass) and 
redshift (e.g. iBlain et a"Dl2003l) : red objects can be either cool 
local galaxies or warmer distant galaxies. What is observed is a 
mis-prediction of the redshift distribution, reflected for example 
in an overestimation of high redshift counts and an underesti- 
mation at later epochs — or vice versa. Several models seem to 
systematically over predict the number of galaxies above z — 1 
in the deep regime. A variety of SED libraries are adopted in 
these model recipes, based on different assumptions and tem- 
plates. These differences in SED shapes and their evolution — 
as well as the implementation of luminosity and density evolu- 



tion for the adopted galaxy populations — indeed produce sig- 
nificantly different results. 

Among all, the Bet hermin et a i1 (l20T0cl) distribution is prob- 
ably the closest to observations, although it presents a signifi- 
cant underestimation of dN/dz at 70 /mi. This model was opti- 
mized taking into account differen tial number c ounts between 15 
Ami and 1.1 mm, in cluding PACS dBerta etalJl2010h and SPIRE 
dOliver et al J 120101) early results, mid-IR luminosity functions 
(LF) up to z = 2, the far-IR local LF, and CIB measurements. 
This success demonstrates the need to include in model tuning 
(ideally by means of a proper automated fit) not only Herschel 
data, but also the detailed redshift information (e.g. evolving LF, 
number counts split in redshift bins, redshift distributions, etc.). 

In summary, while several pre-Herschel backward evolution- 
ary models already provide reasonable descriptions of the to- 
tal counts in at least some of the PACS bands, they tend to fail 
in the synopsis of all counts and in particular redshift distribu- 
tions. The new Herschel dataset is a treasure box, allowing to 
explore the evolution of counts and CIB with unprecedented de- 
tail and much deeper than previous data, on which models were 
calibrated. Significant modifications on model assumptions for 
SEDs and/or evolution will be needed for a satisfactory fit to 
this new quality of data. 

3. Level 2: stacking of 24 //m sources 

Limiting the number counts analysis to individually-detected 
sources only, we miss a significant fraction of the informa- 
tion stored in PACS maps. It is possible to recover part of 
this information by performing stacking of so urces from deeper 
data (typically at shorter wavelengths, e.g. iDole et ail [2006; 
iMarsdenet ail 120091: iBefhermin et alJl2010ah . The 24 /mi cat- 
alog in the G OODS -S field, extending down to ~20 /Jy 
( Mag nelli et al.l2009h . provides the ideal priors to perform stack- 
ing of faint sources on PACS maps. 

The aim here is to transform 24 /mi faint number counts into 
PACS counts, given the average PACS/24/mi colors of galaxies. 
The overall procedure consists in building such colors as a func- 
tion of 24 /mi flux, including both sources detected by PACS and 
un-detected ones. Then the derived colors are used to transform 
24 /mi counts into PACS counts, probing the faint end of the 
distribution, below the PACS detection limit. 

To reach this goal, we select 24 /mi sources not detected 
in PACS maps, and bin them by their 24 /mi flu x. Follo wing 
the standard te c hnique described by IDole et aTl d2006l) and 
Bethermin et alJ d2010al) . in each S(24/mi) flux bin, we pile- 
up postage stamps at the position of these sources and pro- 
duce stacked frames at 70, 100 and 160 /mi. We measure 
the flux density of the stacked sources by performing PSF- 
fitting in the same way a s for the individually-detected ob- 
jects (see Berta et al] |2010l) . Uncertainties on the stacked fluxes 
are computed through a simple bootstrap procedure. Stacking 
was performed on PACS residual images, i.e. after removal of 
individually-detected sources, and the fluxes of the latter were 
then added back to stacking results. Finally, we checked that no 
significant signal was detected when stacking at random posi- 
tions, and that stacks of PACS-detected sources retrieved their 
actual total summed flux. Flux corrections due to losses during 
the high-pass filtering process were tested via simulations and 
tu rned o ut to have a negligible effect on our results (see Lutz et 
al. 201 1 for a description). 

We derived average flux densities simply dividing by the to- 
tal number of sources in each 24 //m flux bin, and then we built 
PACS/24/mi average colors. Figure|6]shows the results. Stacking 
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Fig. 6. Stacking of 24 pm sources on PACS maps in GOODS- 
S. Top panel: PACS/24-yum colors of sources, as measured from 
individual detections (black dots) and detections + stacking (col- 
ored symbols). Black solid lines and grey shaded areas represent 
average colors and dispersion of PACS-detected sources. Solid 
purple lines are polynomial fits to the stacked points. For each 24 
pm flux bin the number of sources detected by PACS and the to- 
tal number of MIPS 24 pm objects ar e quoted. Bottom panel : 
Completeness analysis, based on the iLe Borgne et alj {2009) 
backward-evolution model; the quoted 24 pm fluxes are in [mJy] 
units. 



of mid-IR sources actually allows us to retrieve missed PACS 
fluxes down to faint regimes, thus recovering the actual average 
colors, that would not be possible to derive otherwise. On the 
bright end, average stacked+detected fluxes and individual de- 
tections progressively converge, and become consistent within 
errors when >55% of 24 pm sources in the given flux bin are 
detected by PACS. 

Stacked points were fit with a polynomial function and then 
used to transform the observed 24 pm differential counts into 
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Fig. 7. Differential number counts dN/dS in GOODS-S. 
Resolved counts (black symbols), results of stacking (blue 
crosses) and P(D) results (red solid line and 3<r grey shaded 
area) are shown. 



PACS c ounts, following the recipe suggested by Bethermin et al. 
(l20T0ah : 



dN 



dS 



dN dS 24 
x 



PACS 



dS 



24 



dS 



(1) 



FIR 



where the right terms are computed in the given 24 pm flux 
bin and the left term is computed at the corresponding S fir = 
f(S 24), as defined by the average color functions. The derivative 
dS 24 IdS fir is estimated numerically. 

Although this method has been widely and successfully used 
by several authors, it is worth to note that stacking results are af- 
fected by completeness limitations. In this specifi c case , a flux 
cut in the 24 pm priors (e.g. the Magnelli et al. 120091 catalog 
reaches 20 piy at the 3cr detection threshold) induces incom- 
pleteness in the PACS stacked fluxes, because red objects (faint 
at 24 pm) are not included in the stack. Consequently, stacking 
of 24 pm sources on the PACS maps provides only a lower limit 
to the PACS number counts belo w a given PACS flux. T he bot- 
tom panel of Fig. [6] is based on ILe Borgne et alj J2009) mock 
catalogs, and shows that, for a 24 pm cut of 20 pJy, the PACS 
stacked counts are 90% complete only above ~0.35, ~0.7 and 
~1.3 mJy at 70, 100 and 160 pm, respectively. To this effect, one 
should add the intrinsic properties of the 24 pm parent catalog, 
which reaches ~90% completeness at S (24) ~ 35 pJy. 

FigureQshows the differential number counts dN/dS and in- 
cludes the results of stacking above the 80% stacking complete- 
ness (big blue crosses), in good agreement with the faint-end of 
the resolved counts presented before (open symbols in Figs. Q 
and [TJ. The result of a power-law fit is reported in Tab. [2] at 
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70 and 100 /mi these slopes are slightly flatter, while at 160 /mi 
slightly steeper, than those of resolved counts. 

4. Level 3: P(D) analysis 

It is possible to extend the study of number counts and 
CIB to even deeper flux regimes, beyond the confusion 
limit, exploiting the so-called "pro bability of defl e ction" 
statistics, or P(D) distribution (e.g. IScheuer & Rvld 1 1957 ; 
Condonll974t[F rancesc hini et al .11 1 9 89L [20 1 Ot lOliver et al.ll997 ; 
Patanchon et al.ll2009b iGlenn et al]l2010l) . In a simplified nota- 
tion, this is a representation of the distribution of pixel values 
in a map: its shape and width are mainly driven by three compo- 
nents: number counts, the instrumental spatial response function, 
and the instrumental noise. Once the three are known, it is possi- 
ble to reproduce the observed pixel flux probability distribution, 
or — vice versa — the observed P(D) can be used to constrain 
the underlying dN/dS number counts. 

The first P(D) applications were carr i ed out at ra- 
dio frequencies (e.g. [ Scheuer & Rvld 119571: ICondonl [l974; 
Francesc hini et alj Tl989), but have subsequently been applied 
to different wavelengths, spanning across the whole electromag- 
netic spectrum (e.g. X-ray, Scheuer et al. Il974| Ba rcons et al. 
1994 1 : sub-mm, Hughe setal. [19981 WeiB et al. l2009b Sco tt et al. 
20 id Patanchon et al. l2009t infrared, Oliver et al Il997l) . Since 
the Herschel satellite was launched, a number of analyses have 
made use of P(D) techniques to describe t he properties of real 
maps (IGlenn et al.ll2010t lOliver et al.l l2010). as well as produce 
model predictions to be compared to actual far-I R observations 
dBefhermin et al.ll2010ctlFranceschini et al.ll2010l) . 

We performed the P(D) analysis on GOODS-S PACS maps, 
modeling differential counts dN/dS as a broken power-law with 
three sections and two nodes. To this aim, we developed a new, 
numeric method to estimate the P(D), given the counts model, 
the observed PSF and the observed noise. The position of nodes 
and amplitude of counts at these positions are free parameters. 
The best combination of parameters reproducing the observed 
P(D) is sought by minimizing the difference between model pre- 
diction and data, via a Monte Carlo Markov Chain (MCMC) en- 
gine. Appendix lAl describes the details of this method. 

Figure[8]shows the fit to the observed P(D) in the three PACS 
bands; Fig.[7]and Tab.fTOlreport on the number counts model pro- 
ducing the best fit to the observed P{D). Note that no constraint 
was set to number counts, but only the observed P(D) distribu- 
tion was fit. It is remarkable to note the very good agreement be- 
tween P(D) results and the resolved number counts at the bright 
end, and over the whole flux range covered by the latter. At the 
very faint end, the number count model derived from P{D) tends 
to diverge, because of degeneracies. This is dramatically evident 
at 70 /mi, where the faintest section in dN/dS is totally unde- 
fined. For this reason, we limit the 70 //m result to ~ 0.35 mJy, 
when computing CIB surface brightnesses in the next Sections. 
At 100 and 160 /mi, the P{D) approach allows the knowledge 
of PACS number counts to be extended down to ~0.2 mJy, one 
order of magnitude deeper than individual detections. The re- 
sulting slopes are consistent with resolved and stacked counts, 
over the flux range in common. 

5. Confusion noise 

The extreme depth reached with the P(D) analysis, and the de- 
tailed view of number counts coming from individually-detected 
sources, make worthwhile to take a digression and study the 
properties of confusion noise in PACS maps. 



Following the lDole et al] (120031) formalism, two different ap- 
proaches can be adopted to describe the confusion noise due to 
extragalactic sources, i.e. the effect of fluctuation s due to the 
presence of point sources in the beam (see a lso iFraver et all 
2006; ICondonl [19741 iFranceschini et al] 119891, among others): 
the photometric and the density criteria. 

Given a detection threshold Si m , the former is derived from 
the signal fluctuations due to sources below S n m , while the lat- 
ter consists in deriving the density of sources detected above 
Sn m and limiting the probability that objects are missed because 
blended with neighbors to a small fraction (e.g. 10%). 

Here we apply both criteria to PACS/PEP data and directly 
measure the confusion noise and limits of Herschel deep extra- 
galactic imaging. 

5.1. Photometric confusion 

Photometric confusion is defined by setting the desired signal-to- 
noise ratio q between the limiting flux S ;„„ and the noise 07. de- 
scribing beam-to-beam fluctuations by sources fainter than S i m , 
i.e. q = S u m /cr c . We adopt a value q = 5 for our analysis. 

Generally speaking, the noise measurable on PACS images 
is given by the combination of the instrumental noise 07 (in- 
cluding photon noise, detector noise and data-processing noise) 
and 07. The instrumental n oise was estimated empirically, sim- 
ilarly to Fraver et al] d2006l) . We build partial-depth maps in the 
GOODS-S field, exploiting only a fraction of the available PACS 
observations. By progressively increasing total exposure time, 
we drew the diagrams in Fi g. [9] Sources w ere detected above 
Sii, n in the usual way (see iBerta et al.ll2010i and Lutz et al. in 
prep.) and subtracted prior to measuring the noise on the maps. 
For short exposure times, noise is dominated by instrumental ef- 
fects, and is proportional to f~ . This allows to estimate 07 and 
extrapolate it to deeper regimes. When a long effective exposure 
time is used, the measured total noise cr T deviates from the r - 5 
trend, testifying that confusion noise plays an increasingly sig- 
nificant role. Since both 07 (after source extraction) and 07 have 
nearly Gaussian distributions, the confusion noise can also be 

approximated by 07 = ^cr\ - crj . We iterated between source 

extraction at different depths and noise measurements, until con- 
vergence at q = 5 was reached. 

No deviation from the instrumental 07 oc r 05 is detected at 
70 /mi, at the depth of PEP GOODS-S observations, while we 
derive 07 = 0.27 and 0.92 mJy at 100 yum and 160 /mi, respec- 
tively (with q = 5). Table ITT1 summarizes the derived values. 

5.2. Density criterion for confusion 

The second criterion to determine the influence of confusion on 
extragalactic observations is defined by requiring a minimum 
completeness level in the detection of sources brighter than S u m , 
driven by the fraction of objects missed because blended to their 
(detected) neighbors. Source number counts provide a neat one- 
to-one relationship between source density and flux, so that the 
flux at which the source density confusion (SDC) limit is reached 
is univocally defi ned, given dN/dS . 

Adopting the Laga che et alj d2003l) definition of beam (i.e., 
Q. - 1.14 x 2 FWHM ), and allowing a probability P = 10% 
that sources are t ightly blended and thus cannot be extracted, 
iDole et alj (120031) derived a SDC limit of 16.7 beams/source. 
In what follows we use this value in order to facilitate com- 
parisons with previous works and future exploitations of this 
PACS confusion limit estimate. Nevertheless, it is important to 
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Fig. 8. P(D) distributions in GOODS-S. The black histogram and error bars belong to the observed P(D), while the red solid line 
and grey shaded area represent the best fit and its 3<x confidence interval. The Gaussian noise contribution is depicted with blue 
dashed lines. 



recall that such threshold depends on the steepness of source 
co unts. Following the g eneral treatment of counts and confusion 
bv lFranceschinil ( 19821) . the density of sources with flux S > S /,,„ 
can be expressed as n(> S n m ) (-i-a)q 2 > wnere a is me slope 
of number counts, and q was defined in Sect. 15. II The canonical 
density of 16.7 beams/source is retrieved for a = -1.8, which 
lies within the range of slopes estimated by power-law fitting of 
PEP data (see Tab. |2j, but the SDC limit can significantly vary 
between 8.3 and 25 beams/source for a — —1.5 to -2.0. 

The depth reached in PACS maps, especially in the GOODS 
fields, is such that the high density of detected sources hinders 
the extraction of fainter objects. In Paper I we have already 
treated the case of GOODS-N, showing that the observed cat- 
alog is already hitting the density confusion limit at 160 /mi. 

At 70 //m, despite deep observations, GOODS-S is far from 
being confused: even at the 3<x limit of 1 . 1 mJy, PACS catalogs 
have a density of 32 beams/source including completeness cor- 
rection. The SED shape in the far-IR is such that, while a large 
number of sources is detected at 100 and 160 /mi at the PEP ef- 
fective exposure time, positive ^-correction strongly suppresses 
the number of detections at 70 //m. The typical 70/160 color of 
galaxies is red enough to require a deeper threshold in order to 
detect a comparable number of sources in the two bands, but 
PACS sensitivity at 70 //m is not sufficient to compensate this 
effect. Based on P(D) results, the 16.7 beams/source SDC limit 
would be reached at ~ 0.4 mJy in this band. 

At longer wavelengths, GOODS-N and GOODS-S 
completeness-corrected number counts and P(D) modeling 
agree in indicating that the SDC limit is reached between 1.5 
and 2.0 mJy at 100 /mi and around 8.0 mJy at 160 /mi (Tab.fTTTi. 
At the GOODS-S 3cr levels (1.1 mJy at 100 /mi and 2.4 mJy 
at 160 /mi) the density of sources in our catalogs is 10 and 5 
beams/ source, afte r comp leteness correction. It is worth to note 
that the lDole et al.l d2003l) SDC limit corresponds to a blending- 
completeness of 90%, while sources can be reliably extracted at 
lower completeness also below the 16.7 beams/source threshold. 
Furthermore, photometric completeness further flattens the de- 
tected counts, thus allowing for sporadic very faint objects to be 
extracted. All values quoted above here are nevertheless based 
on completeness-corrected number counts. For comparison, the 
depth reached by P{D) at 100 and 160 /mi digs deep into density 
confusion, down to levels corresponding to 2-4 beams/source. 

6. Discussion: the Cosmic Infrared Background 

Since its discovery by COBE dPuget et al.l Il996t lHauser et all 
1998), several attempts have been carried out to derive the sur- 



face brightness of the CIB, based on direct measurements, inte- 
gration of number counts, statistical analyses, and y — y opacity 
to TeV photons. The PACS data exploited so far can now be 
used to derive new, stringent lower limits to the actual far-IR 
background. 

6.1. CIB measurements: state of the art 

The integral of number counts provides a lower limit to the 
cosmic IR background, to be compared to direct measurements 
from COBE map s. Here we adopt the DIRBE measurements by 
iDole et al] (I2006I) : 14.4 ± 6.3, 12.0 + 6.9, and 12.3 ± 2.5 [nW 
rrT 2 sr _1 ], at 100, 140, and 240 //m, respectively. Performing 
the geometrical aver age between t he up per and lower limits to 
the CIB given by the lFixsen etaf] (1 19981) fit to the FIRES spec- 
trum, one obtains a 160 /mi CIB brightness estim ate of 12.8+6.4 
[nW irT 2 sr _1 ], cons istent with l Dole et al.l (12006) measurements 
at 140 and 240 /mi. IDole et all d2006l) thoroughly describe the 
calibration and uncertainty details for CIB direct measurements. 
The most significant problems arise from the effective bright- 
ness of zodiacal light, currently known only within a factor of 
tw o accuracy. 

Odegardetal. (2007]) used the Wisconsin Ho- Mapper 
(WHAM) Northern Sky Survey as a tracer of the ionized 
medium, in order to study the effects of the foreground interplan- 
etary and Galactic dust on DIRBE CIB measurements. The au- 
th ors find CIB surfa ce brightness values in agreement with those 
by IDole et aljj2 006). onc e reno rmalized to the FIRAS photo- 
metric scale. Ijuvela et al] d2009l) derive the CIB surface bright- 
ness at 170 /mi from ISO maps, obtaining a value ~ 1.5 times 
higher than COBE results, with 30% systematics and 30% statis- 
tical uncertainties. Based on recent observations of the AKARI 
Deep Field South (A DF-S), carried out with this satellite at 65, 
90, 140 and 160 ^m. iMatsuura et al.l d2010l) detected and mea- 
sured the absolute brightness and spatial fluc t uation s of the CIB, 
deriving values consistent to our lDole et al.l (120061) 160 /mi ref- 
erence, but detecting a significant excess at 90 - 140 /mi. 

Finally, the most recent estimate of the far-IR CIB comes 
from the fluctuation analysis carried out by Penin et al. (1201 ll) 
on Spitzer 160 /mi maps of the SWIRE ELAIS-N1 field and re- 
proce ssed IRAS 60 and 100 //m data dMiville-Desch enes et al. 
12002b . The newly-derived CIB brightness is 6.6 + 2.7 and 14.4 + 
2.3 [nW mT 2 sr _1 ], at 100 and 160 /mi, respectively. The 160 
/mi measurement is consistent with the previous COBE interpo- 
lation, but the 100 /mi value is almost halved. 

No direct measurement is available at 70 /mi, except for the 
Miville-Desch enes et al.l(l2002l) fluctuation analysis on the IRAS 
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Fig. 9. Noise in PACS GOODS-S maps at 70, 100, 160 /mi (top, 
middle, bottom panels), as a function of exposure time. The 
dashed line represent the trend 07 oc t 5 followed by instru- 
mental noise; while the dotted line is crj obtained by summing 
07 and the confusion noise <x c in quadrature. Squares mark the 
expected pure-instrumental noise at the exposure times consid- 
ered, while circles denote the measured <x r . 



60 //m map (v/ v = 9.0 [nW irT 2 sr~']) and the AKARI estimate 
at 65 /mi (vl v = 12.4 ± 1.4 ± 9.2 [nW irT 2 sr~'], including statis- 
tical and zodiacal-light uncertainties, Matsuura et al. 2010). 

Further constraints on the actual value of the CIB come from 
the cosmic photon-photo opacity: very high-energy photons suf- 
fer opacity effects by y—y interactions with local radiation back- 
grounds, produ cing particle pairs (e.g. lFranceschini et alj [2008; 
Nikishov 1962). Because of the large photon density of the cos- 



mic microwave background (CMB), any photon with energy 
e > 100 TeV has a very short free path, and extragalactic sources 
are undetectable above this energy. Less energetic photons suf- 
fer from the opacity induced by extragalactic backgrounds, other 
than the CMB, including the CIB. Observations of absorption 
features and cut-offs in the high-energy spectra of BLAZARs 
have been s uccessfully used to pos e upper limits to the intensity 
of the EBL. iMazin & Ra ue (2007), among many others, assem- 
bled a compilation of 1 1 TeV BLAZARs (all those known at the 
time of their analysis) and combined them obtaining constraints 
to the EBL over the whole 0.44-80 /mi wavelength range. 

6.2. The total CIB in PEP 

Integration of number counts was performed over as wide a 
flux range as possible, i.e. using the combined counts pre- 
sented in Sect. [2] including all the four fields analyzed so far: 
GOODS-S, GOODS-N, Lockman Hole and COSMOS. To this 
aim, GOODS-S data were extended down to the 3<x detection 
threshold, thus reaching 1.2 mJy at 70 and 100 //m and 2.0 mJy 
at 160 /mi. At the bright side, COSMOS allows the integration 
to be carried out to 140 mJy at 100 /mi and 360 mJy at 160 /mi. 
Including completeness corrections, the derived CIB locked in 
resolved number counts is vl v = 7.82 + 0.94 and 9.17+0.59 [nW 
irT 2 sr" 1 ], corresponding to 54+7% and 72+5% of the COBE di- 
rect measurements at 100 and 160 /mi, respectively. If the Penin 
et al. d2011l) 100 /mi value were taken as reference, then PACS 
sources would produce ~100% of the CIB. The resolved CIB at 
70 urn is v/ v = 3.61 + 1.12 [nW m~ 2 sr 1 ]. 

The PEP fields are missing the very bright end of num- 
ber counts, that can be probed only over areas even larger than 
COSMOS. This lack of information is particularly relevant at 
70 /mi, because this band was employed only for GOODS- 
S observations. We therefore add previous bright counts mea- 
surements to our data, in order to derive the CIB brightness 
emitted above PEP flux limits, integrated all the way to +00. 
The available data taken into account here are the 70 and 
160 /mi Spitzer counts bv iBethermin et ail d2010al) . built over 
a more than of ~50 deg 2 and ext ending to ~1 Jy, and 100 
/mi ISO and IRAS number cou nts (Rodighiero & France schinil 
2004 : iHeraudeau et"al1 12004 iRowan-Robinson et alJ 
Efstat hiou et al.ll2000UOliver et al.lll992tlBertin et al.l l 
tending to ~60 Jy. Beyond the flux range covered by these past 
surveys, we extrapolate the number counts with an Euclidean 
law, dN /dS oc S ~ 2 ' 5 , normalized so to match the very bright end 
of observed counts, and extended to infinity. 

Figure [10] shows the cumulative CIB surface brightness as 
a function of flux, as derived from the combination of PACS, 
Spitzer, ISO and IRAS data. The curve describing the resolved 
CIB (red circles) rapidly converges to the COBE measurements 
at 100 and 160 //m, over the flux range covered by PACS/PEP. 
Our data extend the knowledge of counts and CIB by over an 
order of magnitude in flux, with respect to previous estimates 
based on individual detections. The depth reached in GOODS-S 
is similar to the one obtained in Abell 2218 through gravitational 
lensing, and the corresponding CIB surface brightnesses are con- 
sistent within the uncertainties. The total CIB emitted above PEP 
3<x flux limits (1.1 mJy at 70 /mi, 1.2 mJy at 100 /mi, 2.0 mJy at 
160 /mi), is vl v = 4.52 + 1.18, 8.35 + 0.95 and 9.49 + 0.59 [nW 
irT 2 sr~'] at 70, 100 and 160 /mi, respectively. Table [T21 lists all 
the values obtained. 

Roughly 10% of the CIB down to the 70 /mi adopted thresh- 
old is emitted by bright galaxies, out of reach for our survey, 
due to the limited volume sampled. On the other hand, only a 
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Fig. 10. Cumulative CIB as a function of flux. Red circles, blue 
crosses and the yellow square belong to completeness-corrected 
counts, stacking and P(D) analysis, respectively. The black line 
solid and grey shaded area are based on power-law fit to the 
resolved number counts. The star symbol indica tes the CIB 
derive d from gravitational lensing in Abell 2218 (lAltieri et al.l 
12011 . The horizontal lines and shades mark the reference di- 
rect measurements of the CIB by iDole et all (120061 at 100 and 
160 nm , including l<x uncertainty) and Miville-Deschenes et al.l 
(2002, at 70 /mi). The green shaded area belong s to previous sur- 
veys carried out with IRAS, ISO and Sp itzer (Oliver et al. 1992; 
Berti n et al.ll997tlE fstathiou et al. 200 0t lRowan- Robinson et a~U 
2004; Rodighie ro & Francesc hini 2004; lH6raudeau et al.l l2004). 



At the very bright end, an Euclidean extrapolation is used (dotted 
lines). 



few percent is produced at fluxes brighter than PEP COSMOS 
upper limit at 100 and 160 /mi. Including 100 /mi ELAIS and 



IRAS data (up to 60 Jy), only a negligible fraction (< 0.05%) is 
missed at the bright end, while beyond Spitzer counts (reaching 
~1 Jy at 70 and 160 //m) the Euclidean extrapolation provides 
roughly 1 -2% of the total CIB . 

The power-law fit to resolved counts (see Sect. [2] and Tab. 
|2]l provides a new estimate of the total expected CIB. We ex- 
trapolate power laws down to 1 .0 //Jy and obtain vl v > 11 .09, 
vl v = 12.61^1 and vl v = 13.63^ | [nW irT 2 sr" 1 ] in the three 
PACS bands. At 70 /mi our data provide only a lower limit, be- 
cause the curve obtained from our power-law fit is not fully con- 
verging at 1 //Jy (see Fig. [TOl , The uncertainty on the 100 //m 
CIB is still large because discordant slopes were found at the 
faint end in the different PEP fields. The results are fully consis- 
tent with COBE data, but PEP and PACS pinpoint the total CIB 
values with unprecedented precision, thanks to the high quality 
of observations, survey strategy, maps, and — last but not least 
— thanks to Herschel capabilities (grey shaded areas in Fig.fTOl 
denote the 3cr uncertainty on power-law fits). 

Stacking results are fully consistent with resolved number 
counts, and the P(D) analysis extends down by another decade 
in flux with the exception of 70 //m, for which we truncated 
the computation of CIB because of a strong divergence in P(D) 
uncertainties (see Sect. @J. Consequently, the lower limit set to 
the CIB by PACS data is further improved: above the flux level 
reached by P{D) statistics, the contributions to the CIB are 4.98, 
9.32 and 11.31 [nW m -2 sr~'] i n the three bands, thus recover- 
ing -65% and 89% of the total IDole et all (120061) CIB at 100 
and 160 //m, respectively. When referred to the total values ob- 
tained through power-law extrapolations, the P(D) analysis re- 
covers 45%, 74% and 83% of the total background at 70, 100 
and 160 //m. 



6.3. CIB contributions from different cosmic epochs 

Thanks to the rich ancillary datasets available in PEP fields, it 
is possible to estimate the amount of CIB emitted at different 
epochs. In order to reach this goal, we integrate the number 
counts split into redshift bins. The same completeness correction 
was applied in each redshift bin, as derived from total counts. 
Table [Dlreports the results for GOODS-S and for the combina- 
tion of GOODS-S, GOODS-N and COSMOS. 

Figure Q~T] shows the fraction of resolved CIB emitted at dif- 
ferent cosmic epochs in the GOODS-S and COSMOS fields. In 
this case, we keep the two fields separate, in order to study the 
details of different flux regimes and aid the comparison to model 
predictions. As stated above, deep observations resolve most of 
the CIB at 100 and 160 //m, and thus give a fairly complete cen- 
sus of the redshift dependence of the background. At the PEP 
3<x detection threshold, more than half of the resolved CIB is 
emitted by objects lying at z < 1. The new results are in line 
with our Paper I analysis, based on shallower GOODS-N detec- 
tions and stacking of 24 /mi sources. As expected, at the depth 
of GOODS-S, galaxies in the highest redshift bin come into play 
and provide a non negligible contribution to the CIB, as high as 
-15% of the resolved amount in GOODS-S. 

For comparison, based on SED template extrapolations from 
15 /mi. lElbaz et al.l (120021) estimated that the contribution of ISO 
15 /zm sources to the 140 //m CIB was of the order of 65%, 
with a median redshift of z = 0.6, and mostly emitted below 
z = 1 . Exploiting Spi t zer 24 //m observations in the COSMOS 
field, iLeFloc'hetall (120091) showed that -50% of the 24 fim 
background intensity originates atz — l. 

The comoving (IR) luminosity (or SFR) density dependence 
on redshift (Madau, Pozzetti & Dickinson 1998, see Gruppioni 
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Fig. 11. Redshift distribution of CIB surface brightness, in GOODS-S (left) and COSMOS (right). Top panels depict the resolved 
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(Mars den et aB2010tlB6thermin et alJl2010cHWiante et alJl2009HLe Borgne et al.ll2009h . See Fig for model lines notation. 



et al. l2010l for a recent determination based on H erschel data) is 
known to peak between redshift z — 1.5 and 3.0. lHarwi 
showed that when transforming the resulting energy generation 
rate into unit redshift interval, the high redshift component is 
strongly suppressed, by a factor (1 + z) 5 ^ 2 in a qo — 0.5 cosmol- 
ogy. Moreover, the effective energy reception rate, thus the back- 
ground observed locally, is further suppressed by a factor (1 + z), 
because bandwidths are reduced. Based on these and other sim- 
ple considerations, it is easy to demonstrate that the distribution 
of the integra ted CIB should be dominated by z < 1 galaxies 
(lHarwitlll999l) . 

The relative fraction among the four bins moves from low 
to high redshift, as wavelength increases: Fig. [12] shows the 
cumulative CIB fraction as a function of redshift, computed 
on the combination of GOODS-S and COSMOS, and includ- 
ing completeness correction. Half of the CIB detected over the 
flux ranges quoted on each panel was emitted below redshift 
z = 0.58, 0.67 and 0.73 at 70, 100, 160 /mi, respectively. 

While at 70 pm ~80% of the resolved CIB was emitted at 
z < 1.0, at 160 pm this fraction decreases to ~55%. This trend 
is mainly due to the positive ^-correction of far-IR SEDs short 
ward of their peak, and will become even more relevant in SPIRE 
and sub-mm bands , where negative ^-co rrection comes into play 
(lOliver et al.ll2010t iMarsden et aDl2009h . Note also that includ- 
ing the contribution from 70 pm unresolved galaxies would shift 
the distribution to slightly higher redshifts (see also discussion 
about 70 pm depth and confusion in Sect. 15.2b . On the other 
hand, adding the bright end of number counts, sampled by the 
observations in the COSMOS field, does not significantly in- 
fluence the high redshift bins in the CIB, because most of the 
contribution goes into filling the lower redshift slices (see Tab. 
[13). The 90th percentiles of the redshift-cumulative CIB fall at 
z = 1.38, 2.03 and 2.20. Excluding the bright end covered by 
COSMOS (only at 100 and 160 pm), we witness only a marginal 
shift of the 90%-light redshift by A(z) < 0.05. 

Performing SED fitting of each individual object de tected in 
GOODS-S in the PACS bands. iRodighiero et all (1201 Oh derived 
IR luminosities (8-1000 pm) of the sources contributing to far- 
IR counts and CIB. Combining the CIB distribution found above 



and the notion of Malmquist bias, it is no surprise that the re- 
solved background is dominated by normal galaxies (L < 10 11 
L ) at low redshift and more luminous sources dominate the high 
redshift bins. Quantitatively, 95% of the resolved 160 pm CIB at 
z < 0.5 is produced by normal galaxies, >90% of the contri- 
bution at 0.5 < z < 1.0 arises from luminous infrared galax- 
ies (LIRGs, 10" < L < 10 12 L Q ) and ultra-luminous infrared 
galaxies (ULIRGs, 10 12 < L < 10 13 L ) provide 50% of the 
background at 1.0 < z < 2.0 and 88% above. Globally, roughly 
50% of the CIB resolved in the three PACS bands in GOODS-S 
was produced by LIRGs. At 160 pm, the remainder is equally 
distributed between normal IR galaxies and ULIRGs, while at 
shorter wavelengths the former prevail. 

COSMOS shallow observations probe the bright end of 
PACS number counts, and hence are limited to lower redshifts, 
but the large area grants a greater detail in the distribution of 
CIB, allowing a much finer binning. At the COSMOS 80% 
completeness limits (9.0 and 20.0 mJy at 100 and 160 pm, re- 
spectively), the resolved CIB peaks at z — 0.3 and exhibits a 
weak secondary "b ump" at z — 0.7 in both bands. Recently, 
iJauzac et al.l (1201 ll) obtained similar results based on stacking 
of 24 pm sources on Spitzer MIPS 70 and 160 pm. These au- 
thors find a dip at z — 0.5 in the differential CIB brightness in 
COSMOS, i.e. the same secondary peak a z - 0.7 found in PACS 
data. This features are certainly due to the z = 0.73 structure 
known in the COSMOS field. Only a small fraction (<10%) of 
the total CIB is locked in z < 0.5 bright objects, which are domi- 
nated by the non-evolutionary component of number c ounts (see 
for example fFranceschini et al.ll2010tlBerta et al.ll2010h . 

The bottom panels in Fig. QT| present the redshift deriva- 
tive of the background surface brightness, and compare it to 
a set of model predictions. The derivative d(vl v )/dz peaks in 
the 0.5 < z < 1.0 bin and then drops at higher redshifts. 
Le Floc'h et al .1 (120091) find a similar behavior at 24 pm, based 
on deep COSMOS data. The decrement is very steep at 70 pm, 
while it becomes milder at longer wavelengths. While over broad 
redshift binning the five models considered are overall consistent 
to the data, at the low redshift bright-end, where detailed infor- 
mation is available, they significantly differ, especially at longer 
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Fig. 12. Cumulative CIB fraction as a function of redshift, as 
obtained summing the contribution of individual sources in 
GOODS-S and COSMOS, and accounting for completeness cor- 
rection. Top, middle and bottom panels belong to the 70, 100, 
160 fim bands, respectively. Red solid lines and dashed blue lines 
mark the 50th and 90th percentiles, i.e. identify the redshifts be- 
low which 50% and 90% of the detected CIB was emitted. Flux 
ranges covered by observations are quoted in each panel. 



wavelengths. Again, this might be a sign of possible differences 
in SED color- and temperature-luminosity assumptions and evo- 
lution. 

The fine details of the comparison between models and data 
in the COSMOS field are affected by the known z ~ 0.7 overden- 
sity. In order to study the impact of these features, the COSMOS 
and GOODS-S catalogs have been cut to a common (overlap- 
ping) flux range. In this case, the d{vl v )/dz derived in the two 
areas are consistent to each other within the errors, although the 
fine structure of GOODS-S cannot be studied. If then the same 
broad redshift binning is adopted, the features due to COSMOS 
large scale structure are washed out, and the trends seen in the 
two fields fully resemble each other. 

Figure [13] summarizes our findings and compares the PEP 
CIB estimates to direct measurements and upper limits from y- 
rays opacity. Our estimate of the background surface brightness 
emitted above the PEP detection and P(D) limits is fully en- 
closed in the COBE l<x error bars. The uncertainties on the CIB 
values derived from PACS resolved counts, P(D), and power- 
law extrapolation are overall significantly smaller than what 
achieved previously by COBE/DIRBE. Far-IR number counts 
are currently in the position to provide a more accurate esti- 
mate of the extragalactic background light than direct photomet- 
ric measurements. 



7. Summary and conclusions 

Using data belonging to the PACS Evolutionary Probe Herschel 
survey (Lutz et al. in prep.), we derived far-IR number counts 
in the GOODS-S, GOODS-N, Lockman Hole and COSMOS 
fields. By employing stacking and P(D) analyses, it was pos- 
sible to exploit the whole information in PACS maps and extend 
the study to very deep fluxes. The wealth of ancillary data in 
the GOODS and COSMOS fields allowed counts and CIB to be 
split into redshift bins. Through integration of number counts, 
we computed an unprecedented estimate of the contribution of 
IR galaxies to the extragalactic background light, across cosmic 
epochs. Among the many analyses carried out here, the main 
results that were presented are: 

- resolved far-IR number counts at 70, 100, 160 fim extend 
from a few mJy to ~200 mJy. In Euclidean-normalized units, 
our combination of four fields tightly defines the bright side 
of counts, the peak at intermediate fluxes, the downturn and 
the slope at the faint-end. Stacking and P(D) analyses push 
the knowledge of counts by a decade fainter in flux, down to 
-0.2 mJy in GOODS-S. 

- based on the observed counts, we derive a photometric con- 
fusion noise estimate of o~ c = 0.27 and 0.92 mJy (obtained 
with q — 5) at 100 and 160 //m respectively, and 16.7 
beams/source density confusion limits of 0.4, ~2.0 and 8.0 
mJy in the three PACS bands. 

- the fraction of CIB emitted at fluxes brighter tha t PEP 3<x 
thresholds is 58% (74%) of the lDole et atl d2006l) reference 
value at 100 yum (160 //m). Using P(D) statistics, this fraction 
increases to ~65% (~89%); this estimate lies well within the 
lcr uncertainty of the direct CIB measurements, with a 3<x 
confidence interval smaller than COBE's by a factor of 2 at 
160 fim. 

- exploiting power-law extrapolation of number counts at the 
faint side, and including past IRAS, ISO and Spitzer results 
and Euclidean extrapolations at the bright end, we derive 
new expectation values for the total CIB at 70, 100 and 160 



fim: vl v > 11.09, v/ v = 12.61 



3$ and vl v = 



13.63 



+3.53 
-0.85 



[nW 



m -2 sr -lj 

- at the current 3<x detection threshold, more than half of the 
resolved CIB was emitted at redshift z < 1 . The half-light 
redshift lies at z = 0.58,0.67 and 0.73 at 70 (GOODS-S), 100 
and 160 fim (GOODS-S and COSMOS combined), respec- 
tively. The balance moves towards higher redshift at longer 
wavelengths: while at 70 fim roughly 80% of the resolved 
CIB is emitted at z < 1.0, at 160 fim the CIB budget is al- 
most equi -partitioned below/above z- 1.0. Roughly 50% of 
the CIB observed at z — and resolved in the three PACS 
bands in GOODS-S is emitted by LIRGs. 

- most of the available evolutionary models fairly reproduce 
PACS total counts, despite the large variety of adopted as- 
sumptions. On the other hand, the detailed high-quality 
PACS data highlight dramatic differences between models 
and severe failures when compared to redshift distributions 
and CIB derivatives. Only models actually tuned taking into 
account redshift information and early Herschel data are suc- 
cessful in reproducing the new PEP dataset, pointing out 
the need to include these constraints in future modeling at- 
tempts. 

Herschel and PACS set a new reference of the brightness of 
the cosmic far-IR background. Expectantly, next analyses and 
missions, rather than trying to recover the missing fraction of 
CIB, envisage the reconstruction of the bolometric background 
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Fig. 13. The cosmic infrared background. Black filled triangles represent the total CIB emitted above the PEP flux limits, based 
on resolved number counts in GOODS-S, GOODS-N, Lockman Hole and COSMOS, evaluated as described in Sect. 16.21 Yellow 
squares belong to the P(D) analysis in GOODS-S. Histograms denote the contribution of different redshift bins to the CIB, ove r 
the flux range covered by GOODS-S. Literature data include: DIRBE m easurements (filled circles, lcr errors, Dole et al. 2006), 
FIRAS spectrum (solid lines above 200 urn, Lagache et al. 1 19991 12000). Fix sen et"a"D (1 19981) modified Black Body (shaded area), 
60 urn IRAS fluctu ation analysis (cross, Miville-Deschenes et al. 2002), and y-ray upper limits (green hatched line below 80 pm 
iMazin & Rauell2007l) . 



from the mid-IR to sub-mm frequencies. The contribution of 
galaxies at different redshift to the CIB is being investigated to 
an increasing details by current surveys. A further study of the 
redshift-dependent CIB, e.g. based on luminosity functions and 
properly including the effects of ^-correction, goes beyond the 
purpose of this paper and is deferred to future works. This infor- 
mation is not only important in the fine tuning of galaxy evolu- 
tion models, but is also of fundamental importance to constrain 
the intrinsic spectra of distant TeV sources, such as BLAZARs, 
whose y-ray photons interact with the CIB. So far, such con- 
straints have been based solely on model predictions, to be con- 
firmed with direct estimates of the CIB as observed from dif- 
ferent epochs. Herschel source statistics have demonstrated to 
provide a precise estimate of the total CIB. Future studies will 
necessitate to improve direct measurements and thus the study of 
foreground contaminants, in order to complete the understanding 
of the total background and pinpoint its measurement. 
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Table 3. PEP 70 fim number counts, normalized to the Euclidean 
slope. 
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Table 1. Main properties of the PEP fields included in the number counts and CIB analysis. 
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Table 2. Power-law fit to PACS differential number counts in the form dN/dS cc S a . Results belonging to resolved number counts 
(i.e. individually detected sources, plus completeness correction) and stacking of 24 /im sources are reported. 
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Table 4. PEP 100 fim number counts, normalized to the Euclidean slope. 
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Table 5. PEP 160 fim number counts, normalized to the Euclidean slope. 
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11.50 








11.50 








7.24 




8.10 


2.36 


1.00 


8.10 


2.36 


225.42 






8.13 








8.13 








4.81 




6.14 


2.43 


1.00 


6.14 


2.43 


283.79 






9.86 








9.86 








6.45 




<3.18 




1.00 






357.27 


























3.14 


2.39 


1.00 


3.14 


2.39 
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Counts from P(D) 



Field 


Flux range 


Slope 


Error 


& band 


[mjy] 


a 


dot 


GOODS-S 70 


0.2-0.7 


-2.33 


+2.13 


v i v / v / iy -j i \j 


0.7-3.1 


-1.70 




GOODS-S 70 


3.1-100 


-2.44 


m 

-0.09 


GOODS-S 100 


0.2-1.1 


-1.23 


+0.02 


GOODS-S 100 
GOODS-S 100 


1.1-5.8 
5.8-100 


-2.09 
-2.25 


-M 

-0.17 


GOODS-S 160 


0.3-5.8 


-1.20 


+0.30 


GOODS-S 160 


5.8-8.5 


-1.28 


m 


GOODS-S 160 


8.5-100 


-2.73 


-0 08 



Table 10. Results of P(D) fit, including 3 <r uncertainties. 



Band 




S Urn 






16.7 beams/src 


70 /im 




0.4 mJy 


100 /mi 


0.27 mJy 


1.5-2.0 mJy 


160 fim 


0.92 mJy 


8.0 mJy 



Table 11. Confusion photometric noise, obtained assuming q = 
5, and 16.7 beams/source density flux limit, for PACS blank ex- 
tragalactic surveys. 



Description 


band 


Flux range 


vl v 


Error 




fim 


[mJy] 


[nW irr 


2 sr 1 ] 


GOODS-S 


70 


1.2-40 


3.61 


±1.12 


GOODS-S 


70 


>1.2 


4.52 


±1.18 


GOODS-S P(D) 


70 


>0.35 


4.98 


+ 16.2 
-2.07 


Power-law 


70 


Total 


>11.09 




All fields 


100 


1.2-140 


7.82 


±0.94 


All fields 


100 


>1.2 


8.35 


±0.95 


GOODS-S P(D) 


100 


>0.2 


9.32 


+5.47 

- + m 

-1 74 


Power-law 


100 


Total 


12.61 


All fields 


160 


2.0-350 


9.17 


±0.59 


All fields 


160 


>2.0 


9.49 


±0.59 


GOODS-S P(D) 


160 


>0.3 


11.31 


+4.00 
-2.43 
+3.53 
-0 85 


Power-law 


160 


Total 


13.63 



Table 12. The cosmic infrared background inferred from PEP 
data. 



Acknowledgements. We wish to thank M. Bethermin, A. Franceschini, 
C. Gruppioni, G. Marsden, S. Niemi, M. Rowan-Robinson, and E. Valiante 
for providing new, extended model outputs for us, and the anonymous ref- 
eree for useful comments. PACS has been developed by a consortium of in- 
stitutes led by MPE (Germany) and including UVIE (Austria); KU Leuven, 
CSL, IMEC (Belgium); CEA, LAM (France); MPIA (Germany); INAF- 
IFSI/OAA/OAP/OAT, LENS, SISSA (Italy); IAC (Spain). This development 
has been supported by the funding agencies BMVIT (Austria), ESA-PRODEX 
(Belgium), CEA/CNES (France), DLR (Germany), ASI/INAF (Italy), and 
CICYT/MCYT (Spain). 



References 

Altieri, B., Berta, S., Lutz, D., et al. 2010, A&A, 518, L17+ 

Aussel, H., Cesarsky, C. J., Elbaz, D., & Starck, J. L. 1999, A&A, 342, 313 

Balestra, I., Mainieri, V., Popesso, P., et al. 2010, A&A, 512, A12+ 

Barcons, X., Raymont, G. B., Warwick, R. S., et al. 1994, MNRAS, 268, 833 

Berta, S., Magnelli, B., Lutz, D., et al. 2010, A&A, 518, L30+ 

Bertin, E., Dennefeld, M., & Moshir, M. 1997, A&A, 323, 685 



Bethermin, M., Dole, H., Beelen, A., & Aussel, H. 2010a, A&A, 512, A78+ 
Bethermin, M., Dole, H., Cousin, M., & Bavouzet, N. 2010b, A&A, 516, A43+ 
Bethermin, M., Dole, H., Lagache, G., Le Borgne, D., & Penin, A. 2010c, 
ArXiv:1010.1150 

Blain, A. W., Barnard, V. E., & Chapman, S. C. 2003, MNRAS, 338, 733 
Brammer, G. B., van Dokkum, P. G., & Coppi, P. 2008, ApJ, 686, 1503 
Bunker, A. J., Stanway, E. R., Ellis, R. S., McMahon, R. G., & McCarthy, P. J. 

2003, MNRAS, 342, L47 
Capak, P., Aussel, H., Ajiki, M., et al. 2007, ApJS, 172, 99 
Chary, R., Casertano, S., Dickinson, M. E., et al. 2004, ApJS, 154, 80 
Condon, J. J. 1974, ApJ, 188, 279 

Cristiani, S., Appenzeller, I., Amouts, S., et al. 2000, A&A, 359, 489 
Croom, S. M., Warren, S. J., & Glazebrook, K. 2001, MNRAS, 328, 150 
Dickinson, M., Stern, D., Giavalisco, M., et al. 2004, ApJ, 600, L99 
Doherty, M., Bunker, A. J., Ellis, R. S., & McCarthy, P. J. 2005, MNRAS, 361, 
525 

Dole, H., Lagache, G., & Puget, J. 2003, ApJ, 585, 617 

Dole, H., Lagache, G., Puget, J., et al. 2006, A&A, 451, 417 

Dole, H., Le Floc'h, E., Perez-Gonzalez, P. G., et al. 2004, ApJS, 154, 87 

Dommguez, A., Primack, J. R., Rosario, D. J., et al. 2011, MNRAS, 410, 2556 

Draine, B. T. & Li, A. 2007, ApJ, 657, 810 

Draper, A. R. & Ballantyne, D. R. 2011, ApJ, 729, 109 

Efstathiou, A., Oliver, S., Rowan-Robinson, M., et al. 2000, MNRAS, 319, 1 169 

Elbaz, D., Cesarsky, C. J., Chanial, P., et al. 2002, A&A, 384, 848 

Elbaz, D., Cesarsky, C. J., Fadda, D., et al. 1999, A&A, 351, L37 

Fixsen, D., Dwek, E., Mather, J., Bennett, C, & Shafer, R. 1998, ApJ, 508, 123 

Fotopoulou, S., et al. 201 1, in preparation 

Franceschini, A. 1982, Ap&SS, 86, 3 

Franceschini, A., Rodighiero, G., & Vaccari, M. 2008, A&A, 487, 837 
Franceschini, A., Rodighiero, G., Vaccari, M., et al. 2010, A&A, 517, A74+ 
Franceschini, A., Toffolatti, L., Danese, L., & de Zotti, G. 1989, ApJ, 344, 35 
Frayer, D. T., Huynh, M. T., Chary, R., et al. 2006, ApJ, 647, L9 
Frayer, D. T., Sanders, D. B., Surace, J. A., et al. 2009, AJ, 138, 1261 
Genzel, R. & Cesarsky, C. J. 2000, ARA&A, 38, 761 
Gilli, R., Cimatti, A., Daddi, E., et al. 2003, ApJ, 592, 721 
Glenn, J., Conley, A., Bethermin, M., et al. 2010, MNRAS, 409, 109 
Grazian, A., Fontana, A., de Santis, C, et al. 2006, A&A, 449, 951 
Griffin, M. J., Abergel, A., Abreu, A., et al. 2010, A&A, 518, L3+ 
Gruppioni, C, Pozzi, E, Andreani, P., et al. 2010, A&A, 518, L27+ 
Gruppioni, C, Pozzi, E, Zamorani, G, Vignali, C, et al. 2011, submitted to 

MNRAS 
Harwit, M. 1999, ApJ, 510, L83 

Hauser, M. G., Arendt, R. G., Kelsall, T., et al. 1998, ApJ, 508, 25 

Hauser, M. G. & Dwek, E. 2001, ARA&A, 39, 249 

Heraudeau, P., Oliver, S., del Burgo, C, et al. 2004, MNRAS, 354, 924 

Hughes, D. H., Serjeant, S., Dunlop, J., et al. 1998, Nature, 394, 241 

Ilbert, O., Capak, P., Salvato, M., et al. 2009, ApJ, 690, 1236 

Jauzac, M., Dole, H., Le Floc'h, E., et al. 201 1, A&A, 525, A52+ 

Juvela, M., Mattila, K., Lemke, D., et al. 2009, A&A, 500, 763 

Lacey, C. G., Baugh, C. M., Frenk, C. S., et al. 2010, MNRAS, 405, 2 

Lagache, G., Abergel, A., Boulanger, E, Desert, & Puget. 1999, A&A, 344, 322 

Lagache, G., Bavouzet, N, Fernandez-Conde, N., et al. 2007, ApJ, 665, L89 

Lagache, G., Dole, H., & Puget, J. 2003, MNRAS, 338, 555 

Lagache, G., Dole, H., Puget, J., et al. 2004, ApJS, 154, 1 12 

Lagache, G., Haffner, L., Reynolds, R., & Tufte, S. 2000, A&A, 354, 247 

Lagache, G., Puget, J., & Dole, H. 2005, ARA&A, 43, 727 

Le Borgne, D., Elbaz, D., Ocvirk, P., & Pichon, C. 2009, A&A, 504, 727 

Le Fevre, O., Vettolani, G., Garilli, B., et al. 2005, A&A, 439, 845 

Le Floc'h, E., Aussel, H., Ilbert, O., et al. 2009, ApJ, 703, 222 

Lilly, S. J., Le Bran, V., Maier, C, et al. 2009, ApJS, 184, 218 

Lutz, D., et al. 2011, submitted to A& A 

Madau, P., Pozzetti, L., & Dickinson, M. 1998, ApJ, 498, 106 

Magliocchetti, M., Santini, P., Rodighiero, G., et al. 2011, ArXiv: 1105.4093 

Magnelli, B., Elbaz, D., Chary, R. R., et al. 2009, A&A, 496, 57 

Marsden, G., Ade, P. A. R., Bock, J. J., et al. 2009, ApJ, 707, 1729 

Marsden, G., Chapin, E. L., Halpern, M., et al. 2010, ArXiv:1010.1176 

Matsuura, S., Shirahata, M., Kawada, M., et al. 2010, ArXiv: 1002.3674 

Matute, I., La Franca, E, Pozzi, E, et al. 2006, A&A, 451, 443 

Mazin, D. & Raue, M. 2007, A&A, 471, 439 

Mignoli, M., Cimatti, A., Zamorani, G., et al. 2005, A&A, 437, 883 

Miville-Deschenes, M., Lagache, G., & Puget, J. 2002, A&A, 393, 749 

Niemi, S., et al. 201 1, in preparation 

Nikishov, A. I. 1962, Sov. Phys. JETP, 14, 393 

Odegard, N, Arendt, R. G., Dwek, E., et al. 2007, ApJ, 667, 11 

Oliver, S. J., Goldschmidt, P., Franceschini, A., et al. 1997, MNRAS, 289, 471 

Oliver, S. J., Rowan-Robinson, M., & Saunders, W. 1992, MNRAS, 256, 15P 

Oliver, S. J., Wang, L., Smith, A. J., et al. 2010, A&A, 518, L21 + 

Papovich, C, Dole, H., Egami, E., et al. 2004, ApJS, 154, 70 



20 



S. Berta et al.: Building the cosmic infrared background brick by brick with Herschel/PEP. 



Zee. 



0.1 

0.2 
0.3 
0.4 
0.5 
0.6 
0.7 
0.8 
0.9 
1.0 
1.1 
1.2 
1.3 
1.4 
1.5 
1.6 
1.7 
1.8 
1.9 
2.0 
2.1 
2.2 
2.3 
2.4 
2.5 
2.6 
2.7 
2.8 
2.9 
3.0 
3.1 
3.2 
3.3 
3.4 
3.5 
3.6 
3.7 
3.8 
3.9 
4.0 
4.1 
4.2 
4.3 
4.4 



COSMOS 
100 yum 160 jim 



1.01e+03 
1.89e+03 
1.59e+03 
2.01e+03 
1.17e+03 
1.07e+03 
8.93e+02 
1.04e+03 
5.99e+02 
7.79e+02 
4.85e+02 
3.36e+02 
3.96e+02 
1.86e+02 
1.98e+02 
1.80e+02 
1.02e+02 
1.50e+02 
1.80e+02 
1.02e+02 
3.60e+01 
1.20e+01 
2.40e+01 
2.40e+01 
3.60e+01 
5.39e+01 

5.99e+00 
1.20e+01 
5.99e+00 
1.20e+01 

1.20e+01 



7.37e+02 
1.33e+03 
9.83e+02 
1.26e+03 
7.91e+02 
7.07e+02 
5.87e+02 
6.23e+02 
3.96e+02 
5.81e+02 
4.32e+02 
2.88e+02 
3.72e+02 
1.92e+02 
1.68e+02 
1.20e+02 
8.99e+01 
1.38e+02 
1.74e+02 
1.26e+02 
7.79e+01 
2.40e+01 
2.40e+01 
3.60e+01 
5.39e+01 
7.19e+01 
2.40e+01 
1.80e+01 
1.20e+01 
1.20e+01 
1.20e+01 
1.20e+01 
1.80e+01 

5.99e+00 



5.99e+00 5.99e+00 



5.99e+00 



GOODS-N 
100 yum 160 fim 



1.82e+03 
1.04e+03 
3.12e+03 
2.34e+03 
2.86e+03 
2.60e+03 
5.20e+02 
3.38e+03 
3.38e+03 
1.82e+03 

7.81e+02 
5.20e+02 
2.60e+02 
2.60e+02 
2.60e+02 



1.04e+03 
2.60e+02 



2.60e+02 



2.60e+02 



1.56e+03 
7.81e+02 
2.86e+03 
1.82e+03 
2.08e+03 
3.12e+03 

1.82e+03 
1.82e+03 
1.82e+03 

5.20e+02 
5.20e+02 
5.20e+02 
7.81e+02 



5.20e+02 2.60e+02 



1.04e+03 
2.60e+02 

2.60e+02 

5.20e+02 
5.20e+02 



2.60e+02 



70 yum 



GOODS-S 
100 yum 



160 yum 



1.54e+03 
2.77e+03 
2.16e+03 
2.77e+03 
3.08e+03 
3.39e+03 
3.70e+03 

9.24e+02 
2.16e+03 
1.23e+03 
3.08e+02 

6.16e+02 

9.24e+02 



3.08e+02 
3.08e+02 



2.77e+03 
3.08e+03 
4.01e+03 
4.31e+03 
3.08e+03 
5.55e+03 
l.lle+04 
1.54e+03 
1.54e+03 
6.47e+03 
2.77e+03 
2.16e+03 
3.08e+02 
6.16e+02 
9.24e+02 
2.47e+03 

3.08e+02 
6.16e+02 
1.23e+03 
1.23e+03 
1.23e+03 
9.24e+02 
9.24e+02 
9.24e+02 
6.16e+02 



3.08e+02 3.08e+02 



2.16e+03 
2.77e+03 
3.39e+03 
2.77e+03 
3.39e+03 
4.62e+03 
8.94e+03 
1.54e+03 
1.54e+03 
5.55e+03 
3.08e+03 
1.54e+03 
6.16e+02 
1.23e+03 
3.08e+02 
2.16e+03 
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6.16e+02 
1.23e+03 
9.24e+02 
1.23e+03 
1.23e+03 
1.23e+03 
1.23e+03 
6.16e+02 

3.08e+02 
9.24e+02 
3.08e+02 



6.16e+02 6.16e+02 3.08e+02 



3.08e+02 
3.08e+02 



Table 6. Redshift derivative dN/dz [deg ~ 2 ] for PACS detected source in COSMOS, GOODS-N and GOODS-S. Catalogs are cut at 
80% photometric completeness, see Fig.|4]for actual flux ranges covered. 
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Appendix A: A cheap P(D) estimator 

In addition to directly detected sources and to stacking of 24 
yum dense catalogs, further information about the shape of PACS 
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Table 7. GOODS-S 70 /vm number counts, normalized to the Euclidean slope, and split in redshift slices. 
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Table 8. PEP 100 jim number counts, normalized to the Euclidean slope, and split in redshift slices. These counts have been obtained 
with a weighted average between GOODS-S, GOODS-N and COSMOS areas. 



number counts comes from the statistical properties of observa- 
tions, probing the counts at even fainter fluxes. The "probability 
of deflection", or P(D) distribution, is basically the distribution 
of pixel values in a map (see Sect. |4j. 



For a large density of sources, the contribution to the P(D) 
is a Poisson distribution with a large mean value, convolved 
with the instrumental noise (typically assumed to be Gaussian). 
Generally speaking, for sufficiently steep number counts at faint 
flux densities, the depth reached with the P(D) analysis is sig- 
nificantly higher than what is provided by individually-detected 
sources, and is a powerful tool to probe counts slopes in cases 
dominated by confusion. 



Given differential number counts dN/ dS , and a beam func- 
tion (e.g. PSF) f(0,(p), describing the point-source spatial- 



response at po sition x, the probability of deflection P(D) can be 



. robability 

written as (e.g. lPatanchon et aDl2009h : 

" (X 



P(D) = T 



exp 



R(x)e ,bJx dx- R(x)dx 



or, alternatively, iso lating 
iFranceschini et al. 2010): 



P(D) 



Jo "'{-J, 



R(x) [1 - cos (2ttcqx)~\ dx\ X 



x cos 



R(x)sin (2nu)x) dx — 2nuD 



dcj 



(A.l) 



the real part only (e.g. 



(A.2) 



where the mean number of source responses of intensity x in the 
beam is given by: 



R(x)= f 
Jo 



dS 



da> 

m~4 



(A3) 
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Table 9. PEP 160 jim number counts, normalized to the Euclidean slope, and split in redshift slices. These counts have been obtained 
with a weighted average between GOODS-S, GOODS-N and COSMOS areas. 
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Table 13. The contribution of galaxies at different epochs to the cosmic infrared background. 



We defer to Patanchon et alj (120091) . iFranceschini et alj 
d2010h . lGlenn et al. I (120101) for a full treatment and description 
of the P(D) formalism. Only in the case of very simple dN/dS 
functional forms and trivial beams (e.g. Gaussian), the above 
equations can be solved analytically. For an effective beam that 
is not strictly positive, or is not azimuthally symmetric, it is nec- 
essary to use the full 2-D beam map. 

For example, the Herschel/PACS beam is characterized by 
3 prominent lobes, which cannot be described analytically. In 
such a case, if one likes to account for the real beam function, a 
numerical treatment of R(x), and possibly P(D) is necessary. 

Finally, in case of real observations, the instrumental noise 
contribution to pixel flux densities must be included. 



A. 1. The numerical P{D) approach 

Given an instrumental PSF f(6, <p) and a number counts model 
dN/dS , it is possible to predict the pixel flux density distribution 
of a map with M pixels, simply rolling random numbers. We 
describe here the principles of the method we developed in the 
frame of the PEP survey. This approach turns out to be of simple 
implementation, avoiding integrations of oscillating functions or 
Fourier Transforms, and it is accessible to any user with a basic 
knowledge of random number generators. Furthermore, it avoids 
any analytic simplification of the beam function, thus allowing 



to employ the real instrumental PSF, regardless of how complex 
its bi-dimensional profile is. 

For a given dN/dS, we describe how a synthetic P(D) is 
computed. The adopted form of number counts and fitting tech- 
nique are described in the next Section. Parameters are varied by 
means of a MCMC engine and for each realization a new P(D) 
is computed and compared to the observed data. 

The input information needed by the P(D) numerical algo- 
rithm are: 

- description of number counts, dN/dS ; 

- beam function, e.g. in form of an observed PSF 2-D map; 

- value of instrumental noise. 

Further secondary inputs, derived from the three primaries 
above, are: 

- the maximum radius r„ MA in pixels from a given source 
within which a pixel is affected by the source itself, e.g. the 
radius at which 99% of the observed PSF flux is enclosed; if 
the pixel scale of the map is a arcsec/pix, the corresponding 
angular radius is ar max ; 

- the total expected number, No, of sources within the solid an- 
gle Q defined by r max , simply given by the integral of number 

counts: N = Q % dS 

- the cumulative probability that a random source has a flux 
lower than S, given by: P(_S) = ^ §; dS' . 
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As mentioned above, the procedure relies on several dice 
rolls, using random numbers to determine the actual number of 
sources expected to contribute to each pixel of the map, the dis- 
tance of the given source to the pixel, the fluxes of sources. 

Step 1. We assume that the distribution of sources in the sky 
is randomly uniform, i.e. we ignore the clustering proper- 
ties of sources in the sky. Clustering would mainly affect the 
P(D) over scales large r than the size of the GOODS-S field 
under exam here (e .g. iGlenn et al ] 120101: IVieroet al.l [2009: 
lLagache et al.ll2007l) . 

Given A^o, the number of sources affecting the given pixel in 
the current realization is a random Poisson deviate for ex- 
pected A^o- This number is recomputed for each pixel and 
each dN /dS realization. 

Step 2. For each source affecting the given pixel, it is necessary 
to determine the fraction /) of flux affecting the pixel. Roll N 
random, uniformly distributed positions, inside r max (defined 
either by r, <p or by x, y coordinates). Based on the position of 
the given source, and the beam function, derive the fraction 
of flux affecting the pixel for each source. 

Step 3. Finally assign a random flux S ,■ to each of the N sources 
affecting the given pixel, according to the input dN/dS . To 
this aim, the probability P(S) (defined above) is used to 
transform a flat random number distribution in the needed 
flux values. 

Assuming that any further background has been subtracted 
from the map, the flux of the given pixel is finally given by 
D — Do + Yn=i fi$h where Do is the term describing pure noise 
(i.e. random flux, in case of no sources at all) and is given by a 
Gaussian deviate as wide as the noise measured on actual PEP 
r.m.s. maps (see Tab.[TJ. 

A.2. Number counts model and minimization 

Num ber counts are mod eled similarly to Patanch on et al.1 {2009) 
and IGlenn et al.l d2010t) . We adopt a simple, parametric model, 
defined by the differential counts dN/dS at a set of flux density 
knots, or nodes. The model is defined within a S m i„ and S max 
lower and upper flux boundaries; outside of this flux range, we 
simply set dN/dS to zero. Between nodes, number counts are 
described as power-laws, connecting at knots. 

We used a Markov-Chain Monte Carlo (MCMC) sampling 
of the parameter space, to explore the likelihood function of the 
model to reproduce the observed P(D). We use a Metropolis- 
Hastings algorithm as our MCMC sampler. In order for this 
sampler to converge efficiently, we draw new steps from univari- 
ate Gaussian distributions, whose widths were given by dummy 
MCMC runs. 

The position of knots and the corresponding amplitude of 
counts are free parameters, and are varied at each MCMC map 
realization. The upper boundary of the covered flux range is 
fixed to 100 mJy, a value driven by the limits of observed, re- 
solved counts in GOODS-S. In order to minimize the number 
of free parameters, and further degeneracies, we limit the de- 
scription of number counts to a broken power-law with three 
sections, for a total of seven free parameters (3 positions and 4 
amplitudes). 

In the past, a standard approach was to subtract bright indi- 
vidually detected sources from the P(D) analysis, and work only 
at the faint side, but it has been shown that a more robust estimate 
of the source counts is obtained if the P(D) is fitted across th e 
whole available flux density range (e.g. lPatanchon et ai1 l2009). 



Therefore, to derive the observed P(D) to be reproduced, we use 
the observed science maps without subtracting any object. 

Section [4] Figs. [8] Q and Tab. [TUl describe the results of the 
whole P(D) analysis. 

A3. Effects of map-making high-pass filtering 

It is necessary to recall that PEP maps have been obtained by ap- 
plying a high-pass filtering process along the data time-line, in 
order to get rid of 1 // noise. This procedure consist in a running- 
box median filter, and induces two different side effects onto 
PACS maps: 1) the background is naturally removed from the 
maps; 2) the PSF profile is eroded and wings are suppressed. 

At very faint fluxes, and high source densities, close to 
and beyond the confusion limit, the contribution of sources to 
the model P(D) basically generates a "background" sheet of 
sources, which shifts the P(D) histogram to a non-zero peak and 
median. In principle, the source component of the P(D) is al- 
ways positive, and the addition of noise generates the negative 
nearly-Gaussian behaviors at the faint side of P{D). The adopted 
high-pass filtering shifts the whole flux pixel distribution func- 
tion to a zero peak/median. In order to take this into account, 
with the current form of PEP maps, we shift the model P(D) to 
the same zero peak, prior of fitting. 

The second consequence of applying an high-pass filter is 
a modification of object profiles. Detected objects are masked 
and excluded from the median filter derivation. Tests were done 
adding simulated sources to the time-lines before masking and 
before high-pass filtering (see Popesso et al. 2011 and Lutz et 
al. 1201 ll) . The result is that the filtering modifies the fluxes of 
masked sources by less than 5%, and those of unmasked point 
sources — below the detection limit — by a factor < 16%. In 
order to understand the impact of this on the P(D) analysis, 
we built simulated maps 20 times deeper than the GOODS-S 
3cr threshold, once using an un-filtered (i.e. masked) PSF at all 
fluxes, and once using an un-filtered PSF above 3<x, but a filtered 
(unmasked) PSF below. The P(D) obtained in the two cases are 
well consistent to each other, with unit median ratio and a <5% 
scatter across all fluxes; no systematic shift, nor skewness are 
observed. 
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